home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / RANDS55.C < prev    next >
C/C++ Source or Header  |  1996-03-26  |  3KB  |  103 lines
  1. /* ============ */
  2. /* rands55.c    */
  3. /* ============ */
  4. #include "defcodes.h"
  5. /* ==================================================== */
  6. /* rands55 - Returns Uniform Random Variate in ( 0, 1 )    */
  7. /* ==================================================== */
  8. #define MBIG    1000000000        /* from Numerical Recipes    */
  9. #undef    MBIG
  10. #define MBIG    (1073741824/2)        /* 2^29                      */
  11. #define MSEED    161803398        /* from Numerical Recipes    */
  12.  
  13. /* The theory behind this 'subtractive' generator can be found in
  14.  * Knuth, Donald E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of
  15.  * The Art of Computer Programming (Reading, Mass.: Addison-Wesley),
  16.  * paragraphs 3.2 & 3.3.
  17.  */
  18.  
  19. #define mod( a, b )    ( a ) % ( b )
  20.  
  21. /* ===================================================================== */
  22. /* Lrand55 - Returns Uniform Random Long Integer Variate in [0,2^29)     */
  23. /* ===================================================================== */
  24. long
  25. Lrand55(restart)
  26. long    *restart;            /* if non-zero, reinitializes    */
  27. {
  28.     static long ma[55];            /* must be 55, see Knuth    */
  29.     static int start = TRUE;
  30.     static int inext, inextp;
  31.  
  32.     int     i, ii;
  33.     long    mj, mk;
  34.  
  35.     if (start || *restart)
  36.     {                    /* initialize last element    */
  37.     mj = mod(MSEED - labs(*restart), MBIG);
  38.  
  39.     ma[54] = mj;
  40.  
  41.     mk = 1;
  42.     /* initialize remainder of table */
  43.     for (i = 1; i < 55; i++)
  44.     {
  45.         ii = mod(21 * i, 55) - 1;
  46.         ma[ii] = mk;
  47.  
  48.         if ((mk = mj - mk) < 0)
  49.         mk += MBIG;
  50.  
  51.         mj = ma[ii];
  52.     }
  53.     /* 'warm up' the generator    */
  54.     for (ii = 0; ii < 3; ii++)
  55.     {
  56.         for (i = 0; i < 55; i++)
  57.         {
  58.         if ((ma[i] -= ma[mod(i + 31, 55)]) < 0)
  59.         {
  60.             ma[i] += MBIG;
  61.         }
  62.         }
  63.     }
  64.     start = FALSE;
  65.     *restart = FALSE;
  66.  
  67.     inext = 0;
  68.     inextp = 31;            /* 31=55-24, 24 & 55 from Knuth    */
  69.     }
  70.  
  71.     if (++inext >= 55)
  72.     {
  73.     inext = 0;
  74.     }
  75.     if (++inextp >= 55)
  76.     {
  77.     inextp = 0;
  78.     }
  79.     if ((ma[inext] -= ma[inextp]) < 0)
  80.     {
  81.     ma[inext] += MBIG;
  82.     }
  83.     return(ma[inext]);
  84. }
  85. /* ================================================================ */
  86. /* Irand55 - Returns Uniform Random Integer Variate in [0,RAND_MAX] */
  87. /* ================================================================ */
  88. int
  89. Irand55(restart)
  90. long    *restart;
  91. {
  92.     return ((int)(Lrand55(restart) >> 14) & RAND_MAX);
  93. }
  94. /* ================================================================== */
  95. /* Drand55 - Returns Uniform Random Double Precision Variate in [0,1) */
  96. /* ================================================================== */
  97. double
  98. Drand55(restart)
  99. long    *restart;            /* if non-zero, reinitializes    */
  100. {
  101.     return(Lrand55(restart) / (double) MBIG);
  102. }
  103.